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establish  the  strong  consistency  of  the  proposed  estimate.  We 
further  propose  a  one  .  step  modification  of  the  non-iterative 
technique.  It  is  observed  in  the  simulation  study  that  the 
proposed  method  performs  better:  than  the  existing  non-iterative 
techniques  for  reasonably  small  sample  sizes.  The  mean  squared 
errors  of  the  proposed  method  reaches  the  Cramer-Rao  lower  bound 
in  many  situations.  We  also  propose  three  different  kinds  of 
confidence  intervals  and  compare  their  performances  by 
simulation. 
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1. INTRODUCTION 

We  consider  the  following  time  series  model ; 

Vt  \  Cos(aj^t)  +  Sin((j^t)  j  +  (1.1) 

Here  y^'s  are  observed  at  equidistant  time  points,  for  t  =  1,  . 

N.  (w^,  .  .  ,  0)^),  (A^,  .  .  ,  A^)  and  (B^,  .  .  ,  B^)  are  the 

unknown  parameters,  are  distinct  real  numbers  lying  in 

(0,7T)  .  M,  the  number  of  signals  is  assumed  to  be  known  apriori. 
(c^)  is  a  sequence  of  real  valued  i.i.d.  random  variables  with 
E(c^)  =  0  and  V{c^)  =  (1.2) 

The  problem  is  to  estimate  the  unknown  parameters  cj  ,  A  and 
k  =  1,  .  .  ,  M  and  cr  .  The  estimation  of  the  parameters 

of  the  model  (1.1)  is  a  fundamental  problem  in  signal  processing 
(Kay  and  Marple;  1981)  and  time  series  analysis.  The  asymptotic 
theory  of  the  least  squares  estimates  (LSE)  for  this  model  has  a 
long  history.  Whittle  (1951,  1953)  obtained  some  of  the  earliest 
results.  More  recent  results  are  by  Hasan  (1982),  Hannan  (1973) 
and  Walker  (1971).  They  formalized  and  extended  Whittle's 
results.  Walker  (1971)  introduced  the  concept  of  an  approximate 
LSE  for  the  model  (l.i).  He  first  estimated  the  frequencies  by 
finding  the  maximum  of  the  periodogram  and  then  computing  the 
estimates  of  the  amplitudes.  The  approximate  LSE  were  shown  to 
be  strongly  consistent  and  the  asymptotic  normality  of  the 
estimates  were  also  obtained.  It  may  be  noted  that  although 
asymptotically  the  approximate  LSE  estimates  are  equivalent  to 
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the  exact  LSE,  for  finite  sample  sizes  the  performance  of  the 
exact  LSE  are  better  than  the  approximate  ones  in  terms  of  lower 
mean  squared  errors  (Kundu  and  Mitra;  1994).  Kundu  (1993a)  was 
the  first  one  to  give  a  direct  proof  of  consistency  of  the  exact 
LSE  for  the  model  (1.1)  under  the  assumption  of  normality  of  the 
error  random  variables,  the  consistency  and  asymptotic  normality 
for  general  error  random  variables  can  be  found  in  (Kundu  and 
Mitra,  1994). 

It  may  be  noted  that  although  the  least  squares  estimates 
are  the  most  desired  estimates,  the  problem  of  finding  the 
estimates  is  well  known  to  be  numerically  difficult.  Rice  and 
Rosenblatt  (1988)  discussed  the  computational  complexities 
involved  to  obtain  the  LSE.  The  model  (1.1)  being  a  nonlinear 
one,  to  obtain  the  LSE  some  sort  of  iterative  search  procedure 
must  be  employed.  Typically,  search  methods  start  from  an 
initial  guess  value  and  then  proceed  by  a  sequence  of 
Gauss-Newton  steps.  For  this  nonlinear  least  squares  problem  it 
turns  out  that  there  are  many  local  minima  with  a  separation  in 
frequency  of  about  N  which  makes  the  stationary  point  to  which 
the  iterative  scheme  converges  extremely  sensitive  to  the 
starting  values.  This  problem  gets  worse  as  the  sample  size 
increases.  It  is  also  observed  (Rice  and  Rosenblatt; 1988)  that 
unless  the  frequency  is  resolved  at  the  first  step  with  order 
o(N  ) ,  the  failure  to  converge  to  the  global  minimum  may  give  a 
very  poor  estimate  of  the  amplitude.  The  problem  becomes 
especially  severe  if  one  is  estimating  the  parameters  of  several 
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harmonic  components  simultaneously,  since  in  that  situation  the 
iteration  is  taking  place  in  a  higher  dimensional  space  with  many 
local  minima.  The  method  of  Walker  (1971)  for  estimating  the 
initial  values  by  finding  the  maximum  of  the  periodogram  turns 
out  to  have  drawbacks.  A  bias  can  arise  for  moderate  sample 
sizes  that  is  appreciable  compared  to  the  standard  deviation 
suggested  by  asymptotic  theory  (Rice  and  Rosenblatt ; 1988) .  The 
initial  values  provided  to  the  search  algorithms  are  thus 
critical.  A  direct  search  of  the  periodogram  at  a  fine  grid  of 
points  substantially  finer  than  that  given  by  the  frequencies 
27ri/N  used  by  fast  Fourier  transform  is  appealing,  but 
unfortunately  has  its  drawbacks  as  well.  Due  to  the  difficulties 
in  obtaining  the  least  scjuares  estimates  several  non-iterative 
techniques  have  been  proposed  in  the  recent  past.  Among  the 
non-iterative  methods  for  estimating  the  frequencies  of  the  model 
(1.1)  in  the  undamped  exponential  model  form,  the  best  known  is 
the  TLS-ESPRIT  method  of  Roy  and  Kailath  (1989) .  Recently, 
another  non-iterative  technique  have  been  proposed  by  Quinn 
(1994) ,  this  involves  the  computation  of  three  Fourier 
components;  the  Fourier  component  at  the  maximiser  of  the 
periodogram  and  at  the  two  adjacent  Fourier  frequencies. 

In  this  article,  we  propose  a  new  non-iterative  method  for 
estimating  the  frequencies  of  the  model  (1.1).  The  proposed 
method  provides  better  frequency  estimates  (in  terms  of  lower 
mean  squared  errors)  than  the  existing  non-iterative  techniques 
for  reasonably  small  sample  sizes.  The  mean  squared  errors  of 
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the  proposed  method  reaches  the  Cramer-Rao  lower  bound  in  many 
situations  in  the  cases  considered  in  our  simulation  study.  The 
proposed  estimate  can  also  be  used  as  an  efficient  initial 
value  for  any  optimization  algorithm  to  obtain  the  least  squares 
estimates.  First  we  transform  the  model  (1.1)  to  an  undamped 
superimposed  exponential  signals  model  and  then  use  extended 
order  modelling  and  singular  value  decomposition  technique  to 
estimate  the  frequencies.  We  call  the  new  estimates,  the  Noise 
Space  Decomposition  (NSD)  estimates.  Once  the  frequencies  are 
estimated  by  the  NSD  method,  the  linear  parameters  can  then  be 
obtained  using  separable  regression  technique  of  Richards  (1961). 
The  proposed  method  is  shown  to  give  strongly  consistent 
estimates.  Since  the  proposed  method  is  strongly  consistent,  a 
further  one  step  modified  estimate  is  also  proposed  which  already 
have  the  same  asymptotic  properties  as  the  exact  LSE.  Some 
confidence  intervals  of  the  frequencies  are  also  proposed. 

The  organization  of  the  paper  is  as  follows;  in  Section  2  we 
describe  the  estimation  procedure,  consistency  results  are 
provided  in  Section  3.  Modified  estimates  are  proposed  in  Section 
4,  different  confidence  intervals  are  discussed  in  Section  5. 
Some  Monte  Carlo  simulations  study  is  presented  in  Section  6  and 
finally  we  draw  conclusions  in  section  7. 


2.  ESTIMATION  PROCEDURE 

Observe  that  the  model  (1.1)  can  be  written  as  a  linear 
combination  of  2M  complex  exponential  terms  in  the  following  way; 
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=  E  exp(i  u^t)  +  E  exp(-i  (j^t)  +  c^,  t  =  1, 


where  C  =  •=•  A  ~  i  B 

r  Z  r  r 


]  -  5( 


(2.1) 


A  +  i  B 

r  r 


and  i  = 


It  is  well  known  (Prony;  1795,  see  Kundu;  1993b  also)  that 
in  the  noiseless  case  there  exists  an  unique  vector  c  =(c^/  •  •  » 


®2M*i  ^  that ; 


E  c  y  =0 

^  k  ■'t+k 
ksi 

for  all  t=0,  .  .  ,  N-2M-1,  IICII  =  1  and  c  >  0. 


(2.2) 


The  unknown  constants  (c  ,  .  .  ,  c  )  are  such  that  the  roots 

1  2M+1 


of  the  polynomial  equation 


c  +  c 
1  2 


+  c  =0 

2M+1 


(2.3) 


are  of  the  form  =  exp(±  iw,^)  •  Thus  if  we  can  estimate  C,  we 
can  estimate  the  unknown  frequencies  oj^'s  using  (2.3).  Now 
observe  that  (Kahn  et.  al;  1993),  the  condition  that  the  roots  be 
purely  imaginary  means  (2.3)  must  factorize  in  the  form 


”  r  2  2  1 

’’o  ^  -  (2-T7px  +  1  J 


(2.4) 


This  implies  that  c  =  c^^  ^  k  =  1,  .  .  ,  2M+1  and  the  roots 

are  A.  =  exp(±  iu)  )  ,  where  2  -  =  2  Cos((j  ). 

Consider  the  following  N-L  x  L+1  data  matrix 


(2.5) 


for  any  positive  integer  L  such  that  2M  ^  L  s  N-2M. 


Let's  denote  by  T  the  L+1  x  L+1  matrix  given  by 
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1  * 

T  =  —  A  A 

A  A 


(2.6) 


where  denotes  the  conjugate  transpose  of  a  matrix  or  a 
vector.  Observe  that  in  the  noiseless  case  the  matrix  T  has  rank 
2M.  Let  the  singular  value  decomposition  of  T  (see  Rao;  pp.  42, 


1973)  be  as  follows; 


A  A^ 

T  =  r  U  U 
itt  ‘  ‘  ‘ 


(2.7) 


''2  '^2  ''2 

where  >  cr^  >  .  .  >  are  the  ordered  eigenvalues  of  T  and 

IS  the  normalized  eigenvector  corresponding  to  (t\  The 
subspace  generated  by  |  U^,  .  .  ,  is  denoted  by  S  and 


that  of  •{  U  , 
1  2M+1  ' 


U  I  i 

L*lj 


is  denoted  by  IN.  We  call  S  the 


signal  subspace  and  IN  the  noise  subspace.  Let  be  any  basis  of 


the  noise  subspace  IN.  We  write 


= 


1.L+1-2M 


b  ...  b 

Ltl.l  L+l,L+l-2M 


(2.8) 


Observe  that  because  of  (2.2),  in  the  noiseless  situation 
there  exists  an  unique  basis  of  IN  which  has  the  following  form 

c  0  ...  0 

^  c  .  .  .  0 


c  c 

M^l  M 
C  C 

M  H* 


(2.9) 
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Now  observe  that  *  '  '  ^l+i]  ^  basis  of 

the  estimated  noise  space.  Our  main  aim  is  to  obtain  a  basis  of 
IN  which  has  the  form  similar  to  (2.9)  and  to  estimate  C  from 
these. 


Let's  partition  the  matrix  as  follows; 


=[ 


L  ♦  1  -  2M  X  L+1 


11  *  ®12  •  ®13 

L+1-2MXIC  L  +  1-2MX2M+1  L  ♦  1  -  2  M  x  L  -  IC-2M 


(2.10) 


for  k=0,l,  .  .  ,  L-2M.  Now  consider  the  matrix 


[ 


L+1-2MXK  L  ♦  1  -  2  M  X  L -K-2M 


Since  the  above  is  a  random  matrix,  it  is  of  rank  L-2M. 
Therefore  we  can  conclude  that  there  exists  an  unique  L+1-2M 
vector  X  *  0,  such  that 

K*1 


KxL  *  I  -  2H 


L-K-2MX  L  +  1-2M 


(2.11) 


Consider  the  2M+1  vector  where 

^  '  ^IC+1.2M+1^  "  ®12 


(2.12) 


By  properly  normalizing  we  can  make  c  >0  and  iic''**il  =  1  for 

,  1 

k  =  0,  .  .  ,  L-2M.  Therefore  we  can  conclude  that  there  exist 


vectors  X  ,  .  .  ,  X,  ,  such  that 

1  L'^  1  *2M 
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I 

I 

I 


Bl  [='l- 


L-2M+1 


]= 


where  c 


k.l 


>  0  and 


IC  II  = 


1 . 1 


A 

c 

0 

0 


1 , 2M+1 


1  for  k=l, 


in  the  noiseless  situation 


0 

A 

C_ 

A 

C 

0 

0 


2,1 


2,2M+1 


L-2M+1,1 


(2.13) 


L-2M+1,2M+1 

L+1-2M.  Observe  that 


C  =  C 

Let  J  be  the  L+1  x  L+1 


2  ;:l-2h+i  ^ 

=  .  .  =  c  =  c 

exchange  matrix  given  by 


J 


0  0  .  .  .  1 

0  0.  .10 

•  ••••• 
1  0  .  .  .  0 


(2.14) 


(2.15) 


Consider  the  matrix  T  given  by 

f  =  J  T  J 


(2.16) 


Observe  that  the  eigenvalues  of  T  and  f  are  same  and  if  x  is  an 
eigenvector  of  T  corresponding  to  the  eigenvalue  X,  i.e. 

Tx  =  Xx  =»JTJJx  =  XJx  =»f  (Jx)  =  X  (Jx) 
then  J  X  is  an  eigenvector  tor  T  corresponding  to  X.  Let^s 


denote  by  fi  the  subspace  generated  by  |  JU 


2M'H 


JU. 


L*1 


}- 


we 


call  the  noise  space  of  T. 

It  can  be  easily  seen  that  in  the  noiseless  situation  there 
exists  an  unique  basis  of  the  noise  space  of  f  of  the  form 


9 


o  > 


0 


0 


= 


(2.17) 


c  c 

1  M 

•*  ; 

c  ^ 

0 


In  this  case  also  our  aim  is  to  obtain  the  basis  of  the 


estimated  noise  space  W,  i.e. 


=  ’ 


:  JU  1 

L*1  J 


(2.18) 


to  the  form  similar  to  (2.17).  Proceeding  exactly  as  in  the  case 
of  IN,  we  reduce  the  basis  to  the  following  form 


A* 

c 

2,1  . 


2,2M+1 

_  ^1,2M+1  ^ 

A  ^  A  ^ 

such  that  for  each  C  =  (c  ,  . 

K  ,  1 


L-2M+1  ,1 


L-2M+1  ,2M+1 


(2.19) 


c  )  ,  k=l , 

k,2M+l' 


,  L-2M+1; 


>0  and  lie  II  =  1.  As  in  the  case  of  IN,  in  the  noiseless 


situation 


A* 

c  =  c  = 

1  2 


=  c  =  c 

L'fl-2M 


(2.20) 


It  is  further  observed  that  G  =  J  ^1' • '  •*L-2M+1  '  ^*®* 
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•  c 

L-2M+1  ,1 


C  ^  * 

2,1.  C 

L-2M+1 ,2M+1 


A* 

C 

1  ,21(+1 


A* 

,  2  M  +  1 

0  .  .  0 


C  0.0 

:  S.i  •  0 


1 , 2M+1  *•  C 

r,  C  .  L-2M+1,1 

0  2.2M+1. 


0  ‘  C 


L“2M+1,2M  +  1 


(2.21) 


Now  observe  that  since  (2.14)  and  (2.20)  are  true,  it  is 
quite  natural  that  any  one  of  the  c"'  for  k=l,  .  .  ,  L-2M+1  or  C* 


for  k  =  1, 


,  L-2M+1  can  be  used  to  estimate  the 


frequencies.  In  fact  the  use  of  C**=  ^(C*^  +  c*)  ;  k  =  1,  .  . 

,L-2M+1  always  ensure  that  the  estimated  coefficients  of  the 
polynomial  prediction  equation  (2.3)  satisfy  the  symmetry 


constraint  and  roots  of 


;!**  „2M  .  „2M-l  , 

c  Z  +  c  Z  +  .  .  +  c  =  0 

If  .  1  k  .  2  k  ,  2M+1 


(2.22) 


are  of  the  form  exp(±itj^)  ,  for  k  =  l,  .  .  ,  L-2M+1.  We  use  all 

^  it  if 

;  k  =  1,  .  .  ,  L-2M+1  to  estimate  cj,  We  take  the  average  of 
all  's  and  use  (2.3)  to  get  the  final  estimate  u  of  u.  We 
call  the  resulting  estimate  the  Noise  Space  Decomposition 

(NSD)  estimates. 
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3.  CONSISTENCY  RESULTS 

In  this  section  we  prove  the  following  result, 

Theorem  1:  Under  the  assumptions  of  the  model  (1.1),  the  estimate 
u  of  (j  obtained  by  the  method  described  in  section  2  is  strongly 


consistent,  i.e. 


To  prove  theorem  1,  we  need  the  following  lemmas 


(3.1) 


Lemma  1:  Let  P  =  ((P,^))  and  Q  =  ( (qj^) )  are  two  Hermitian  m  x  m 
matrices  with  spectral  decompositions 


p  =  r  5  p  p 

1=1 

a  -  rx,  q,q* 

1=1 


5  fe  5  s 
1  2 


X  a  X  a 
1  2 


(3.2) 


where  S^'s  and  are  eigen  values  of  P  and  Q  respectively  , 

and  are  orthogonal  normalized  eigenvalues  associated  with  5^ 
and  Xj  respectively  for  i=l,...,m.  Further  assume  that 


=  X  =  A,  ,  n  =0<n  <  .  .  <n  =  p  ; 
n  n  o  1  s 

h 


A,i^>  ••  ^  A 

12  S 


and  that 


i,lC”  1,  *.  ,  m 


then  there  is  a  constant  M  independent  of  a  such  that 


(1)  < 


1  “  1  ^  •  •  /  in 


(2)  E  p  p  =  Z  q  q  + 

i=»n  i=n 

h-l  h-l 


(3.3) 


with  c<*'>=  ((c(jj>)), 


|c<jj>  I  s  Ma 
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Proof;  It  follows  from  Von-Neumann's  (1937)  inequality,  see  also 
Bai,  Miao  and  Rao  (1991). 

Lemma  2: 


i  A*A  =  <7^  +  n"->  D  n“-'*  + 


o[/I^ 


logN 


a.  s 


(3.4) 


where  D  .  diag  {  |c^| . |c2|,  |  dJ  j . jD^I  I  and 

2MX2M  "  *  'M'j 


L+1  X  H 


e  1 


-i  (L+l)(j 
e  '  1 


-icj 

e  H 


e  1 


g-i(L+l)(j^  gi(L+l)(a^ 


lU) 

e  M 


Proof  :  We  have 


H  A*A  -  T  =  i  ((t^)) 


with  the  following  renaming  of  the  parameters  of  (2.2)  as 

(  Cj  for  i=l,  .  .  ,  M 


a  = 

i 


i  for  i=M+l,  .  2M 


(J 


-0) 


I-M 


for  i=i,  .  .  ,  M 
for  i=M+l,  .  .  ,2M 


(3.5) 


(3.6) 


(3.7) 


we  have  the  following 

N-L-l 


1  4.  1  „  - 

N  *.k  ■  ff  ^  y 


1  =0 
N-L-l /2M 


1+1  ^1+k 


ilo  S+.J  (j^a^exp(i/3^(l+k))  + 

N-L-l  p  2M  2 

“  N  11^  l“ul  )  +  E  a  exp(i(/3  (l+k)-/3  (1+i))) 

1  «o  ^  u  =  l  u^v  ''  “ 


+  c 


'*4  Z^exp(l@jn.k))  . 
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+  c  ,  e  .  1 

1+1  l+k  J 

=  +  i?2  +  ^^3  +  ^^4  +  (say)  (3.8) 

Observe  that  =  0  (1)  and  i?3  =  o[  /l£0^  j  a.s  and 
^4  =  o[  /  ^Pg^^Qgl  j  a.s  (see  Petrove  (1975,  page  375)  and  by 
the  law  of  iterative  logarithm  of  M  dependent  sequence 


^5  = 


f°( 


/  log  logN 
N 


a^+  0 


this  proves  lemma  2 . 


(/^ 


logN 


if  is^k 
if  i=k 


Lemma  3.  Let  gj^(x)  be  a  sequence  of  polynomials  of  degree  k,  with 
roots  x*”’ ,  .... 


x^"’  for  each  n.  Let  g(x)  be  a  polynomial  of 


degree  Q  ,  with  roots  x^,  . .  x^,  Q  s  k.  If  gj^(x)  - >  g(x) 

as  n  - *  00  then  with  proper  rearrangement  the  roots  of  g  (x)  , 

converge  to  the  roots  of  g(x) ,  ie  to  x^. 

Proof  :  See  Bai  et  al  (1986) 

Proof  of  Theorem  1:  From  Lemma  2  it  follows  that; 


T 

0-2- 

pv 

n"-’*  =  s 

J.  T 

L+1 

D 

and 

T  -  *•“  ) 

<T^I  +  J 

U+1 

D 

J  =  s* 

J 

Observe 

that  the  eigen 

values  of 

S  are 

of  the  form 

\  fcX  a  . 

&  A  ^  Ik 

_  _ 

1  2 

2M  2M+1 

• 

•  •  ^ 

L+l 

since  0 

(L)  _ 

DO  IS 

of  rank 

2M. 

Let  the 

(3.9) 


decomposition  of  S  be 

‘-+1  * 
S  =  A,  s. 


(3.10) 


where  Sj^  is  the  normalized  eigenvector  corresponding  to  the 
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eigenvalue  and  s^'s  are  orthogonal  to  each  other  . 
using  lemma  1 


Therefore 


I  »i  «i 

l  =  2M+l 


a,s  * 

- '  ^  ®i  ®i 

1=2M*1 


(3.11) 


(3 . 10)  implies  that  the  vector  space  generated  by  ■|^2m+1'  *  l|^ 

converges  to  the  vector  space  generated  by  /  .  ..s  I  . 

^  \  2M+1'  '  L+lj 

Now  the  former  one  has  a  unique  basis  of  the  form 


1 . 1 


0 

A 

c 


2.1 


1 , 2M+1  ■*> 

0  2,2M+1 


0 

0 


0 

0 


L-2M+1,  1 


L-2M+1,2M  +  1 


(3.12) 


with  c 

k,i 

k  =  1,  . 

the  form 


>  0  and  II  c’^ll  =  1  where  =  (  c 


k,l'  *  *  '  ^k,2M+i^ 


L-2M+1,  whereas  the  later  one  has  a  unique  basis  of 


c 

1 

0 

0 

c 

C 

0 

2 

1 

C 

c 

.  2 

0 

0 

C 

2M^1 

c 

1 

0 

0 


0 

0 


2H4'1 


(3.13) 


Where  c^  >  0  and  ii  C 

This  implies  that 


=  1  with  c  =  (c  .  .  .  . ,  c  ^  )  . 

1  2M+l' 
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c 


for  k 


1,  .  .  L-2M+1 


(3.14) 


Similar  analysis  for  T  shows  that 

>  c  for  k  =  1,  .  .  L-2M+1 

Thus  we  have 

^  >  C  for  k  =  1,  .  .  L-2M+1 


(3.15) 

(3.16) 


Therefore  from  Lemma  3  we  can  say  that  the  roots  obtained  usin^ 

A** 

Cj^  are  consistent  estimators  of  u  's  for  all  k  =  1,  .  L-2M+1. 


4.  MODIFIED  NOISE  SPACE  DECOMPOSITION  METHOD 
It  may  be  noted  that  the  model  (1.1)  is  a  nonlinear 
regression  model.  So  the  least  squares  method  should  be  the 
preferred  estimation  technique  at  least  as  far  as  efficiency  is 
concerned.  It  is  observed  that  for  this  particular  nonlinear 
least  squares  problem  there  are  many  local  minima  with  a 
separation  in  frequency  about  N“‘  which  makes  the  point  of 
convergence  of  iterative  scheme  extremely  sensitive  to  the 
initial  values.  The  need  for  non— iterative  estimation  procedures 
thus  arise  mostly  from  the  fact  that  good  starting  values  are 
needed  in  numerical  optimization  of  the  residual  sum  of  squares. 
In  some  cases,  for  instance  when  the  signals  have  to  be  detected 
on  line,  the  iterative  least  squares  method  might  be  too  time 
consuming.  If  the  motivation  for  the  use  of  a  non— iterative 
technique  is  the  latter,  the  following  improvement  of  the  NSD 
estimates  can  be  suggested. 

It  is  well  known  (Harvey;  1981,  Ch.  4.5)  that  when  a  regular 


(differentiable  upto  third  order)  is  maximized  through 
the  Newton-Raphson,  scoring  or  a  related  algorithm,  the  estimates 
obtained  after  one  single  round  of  iteration  already  have  the 
same  asymptotic  properties  as  the  exact  least  squares  estimates. 
This  holds,  if  the  starting  values  have  been  chosen  /TT  - 
consistently.  Now,  since  the  NSD  method  is  strongly  consistent 
we  combine  the  NSD  method  with  one  single  round  of  scoring 
algorithm.  This  way  the  asymptotic  error  variances  should  (in 
theory,  at  least)  coincide  with  the  asymptotic  variance 
iJ^i^®spGctive  of  the  distributional  form  of  the  error  term.  We 
C3ll  the  resulting  estimates  obtained  after  one  round  of 
iteration  with  NSD  as  starting  values,  the  Modified  Noise  Space 
Decomposition  (MNSD)  estimates. 

One  way  of  implementing  this  idea  would  be  the  following; 

Let  us  write  the  model  (l.i)  in  the  vector  form 

y  =  A(q)  a  +  c  (4.1) 

where  A(tj)  =  Aj^((j)  A2(u)  .  .  .  Aj^((j)  j  with  Aj^((j)  =  |  Cos(u^k) 

Sin(Wjk)  ....  Cos(£j^k)  Sin((j^k)  j  ,  a  =  ^  j , 

=  =  K  •  •  ■  s)'- 

Now  consider  the  concentrated  residual  sum  of  squares 
Y  [l  -  P^(w)j  Y  =  y*  -  A((J)  Ja*((j)A((j)  j  A*((J)J  Y  (4.2) 

To  obtain  the  least  squares  estimates  first  (4.2)  can  be 
minimized  with  respect  to  ta  and  then  the  estimate  of  a  can  be 
obtained  using  linear  regression  technique.  For  details  see 
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Kundu  (,1993b)  .  We  obtain  the  MNSD  estimates  after  one  step 
minimization  of  (4,2)  using  NSD  estimates  as  starting  values. 


5.  CONFIDENCE  INTERVALS 

In  this  section  we  propose  different  confidence  intervals 
for  the  frequencies.  We  propose  an  asymptotic  confidence 
interval  and  two  bootstrap  confidence  intervals. 


5a.  Asymptotic  Confidence  intervals 
In  this  subsection  we  discuss  the  confidence  intervals  for 
the  frequencies  based  on  their  asymptotic  distribution.  It  may 
be  observed  that  (Kundu  and  Mitra;  1994)  that  the  asymptotic 
distribution  of  the  exact  LSE  of  the  frequencies  is  of  the 
following  form 


N 


-3/2 


r  w  -  w  1 

■  wf  0, 

24  cr^  '1 

L  k  k  J 

\  mm  1  ^ 

(A^B^  J 

distribution 


(5a. 1) 
of  the 


which  eventually  coincides  with  the 
approximate  LSE  proposed  by  Walker  (1971). 

Based  on  (5a.l),  the  following  100(l-a)%  confidence  interval 
for  is  proposed  by 


r " 

\  (ji  -  X  , 

(M 

<  b 

CM 

1'"  f 

24  cr  ] 

1/2  1 

J  (5a. 2) 

L  k  a/2 

^  N^(A^+B^) 

k  k ' 

^(A^+B^)  J 

k  k ' 

Since 

the  MNSD  estimates  proposed 

in  section 

4  has  the  same 

asymptotic 

properties 

A 

as  the  exact 

LSE,  we 

take  the  MNSD 

5b,  Percentile  Bootstrap  Confidence  intervals 
In  this  subsection  we  construct  the  percentile  bootstrap 
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intervals  for  using  the  method  suggested  by  Efron 

(1982) . 

Suppose  we  have  a  sample  of  size  N;  y  ,  .  .  y  coming 

1  M 

from  (1.1).  We  propose  the  following  algorithm  to  obtain  the 
confidence  intervals 

(1)  Estimate  (o)^,  .  .  (j^)  from  y^,  .  .  .,  y^  using  MNSD 

method. 

(2)  Estimate  y,-y,,  i=l,  .  .  ,  N,  where 

y,  =  E  A  Cos((jt)  +  i  Sin((jt)  1. 

(3)  Draw  a  random  sample  of  size  N  from  |  .  .  .  ,  j.  with 

replacement,  let  it  be  |  c^,  .  .  .  ,  c  I  . 

(4)  Obtain  bootstrap  sample  y*,  .  .  .,  y*;  where 

y,  ~  '  i^'i^  .  .  ,  N. 

1 

(5)  Estimate  (u^,  .  .  (j^)  from  y^ ,  .  .  . ,  y*  using  MNSD 

method.  Denote  it  by  u*  ,  K=l,  .  . ,  M. 

(6)  Repeat  the  steps  (3)  to  (5)  NBOOT  times. 

Order  these  NBOOT  estimates  corresponding  to  each 
Estimate  (a/2)  by  NBOOT  a/2th  order  statistics  and 

A 

NBOOT  (l-a/2)th  order  statistics  for  each  set  of 
cj  ^  and  claim  that  (Lp^(a/2),  ^p^(oc/2))  to  be  the  100(l-a)% 
percentile  bootstrap  confidence  interval  for  u  . 

k 

5c.  Bootstrap“t  Confidence  intervals 
In  this  subsection  we  construct  the  bootstrap-t  confidence 
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intervals  based  on  the  method  suggested  by  Hall  (1988) . 

We  propose  the  following  algorithm  for  computing  the 
bootstrap“t  confidence  intervals, 

(1)  From  the  original  sample  y^,  .  .  ,  y^  estimate  (w^,  .  . 
by  MNSD  method  and  then  (A^,  .  .  ,  A^)  and  (B^,  . 
by  separable  regression  technique  and  estimate  cr^  by 

A2  1  ^  (  A  '  2 


I  =l 


I  1 


(2)  Estimate  the  errors  by  c^=  yj“yj/  i=l,  .  . 

(3)  Draw  a  random  sample  of  size  N  from  | 

{  Cb,  .  .  .  ,  Cb  I 

'•1  N  '' 


N. 


} 


with 


replacement,  let  it  be 


(4)  Obtain  bootstrap  sample  y*, 


w  y„ ;  where 

N 


=  y,  +  s  .  N- 

I 

(5)  Estimate  (tj  ,  .  .  w  )  from  y*  .  .  .  ,  y*  using  MNSD 

"1  fl 

method,  denote  it  by  u  and  also  the  estimate  of  as  cr^. 

*  B 

(6)  Obtain  for  each  i=l,  .  M 


T 

i 


, —  ^  ^ 
vN  (u  -cj  ) 

_ ^  k 

A 

a 


B 

(7)  Repeat  the  steps  (3)  to  (6)  NBOOT  times. 

(8)  For  each  order  the  NBOOT  number  of  T  's.  Estimate 

L^g(tt/2)  by  <j^+  vN  cr  I  NBOOT  a/2  th  order  statistics  from 


Tj's  j  and  U^b(“/2)  by  /IT  a  |  NBOOT  (l-a/2)  th  order 

statistics  from  T^'s  j.  Now  claim  that  (i^^ia/2)  ,  lJ^^(,a/2)) 
to  be  the  100(l-a)%  bootstrap— t  confidence  interval  for 
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6,  MONTE  CARLO  SIMULATIONS 

We  have  perfonned  Monte  Carlo  simulation  study  to  ascertain 

the  behavior  of  NSD  and  MNSD  estimates  and  to  compare  it's 

performance  with  the  best  known  non— iterative  techniques  for 

moderate  sample  sizes  and  different  ranqes  of  the  error  variances 
2 

^  •  All  these  simulations  have  been  done  on  the  HP— 9000  computer 
at  the  Indian  Institute  of  Technology  Kanpur,  using  the  IMSL 
random  deviate  generator. 

We  consider  the  following  model  of  one  harmonic  component; 

y^  =  1.5  Cos((jt)  +  1.5  Sin((jt)  +  ,  t=l,  25  (6.1) 

The  error  random  variable  (c^)  is  white  and  Gaussian  with 
variances  .  The  frequency  u  is  taken  to  be  0.257T,  O.SOir,  0.75n 
in  three  different  sets  of  simulations.  For  each  u,  100 
independent  trials  using  different  sequences  are  performed. 
The  variance  of  the  error  random  variables  is  varied  from  0.01  to 
1.5.  In  each  case  we  computed  (u^,...,  uj  by  NSD,  MNSD, 
TLS-ESPRIT  and  Quinn's  (1994)  method.  For  each  u,  we  computed 
the  average  estimates  and  the  mean  squared  errors  (MSE)  over  100 

replications  and  also  the  corresponding  Cramer-Rao  lower  bound 
(CRLB) . 

It  is  observed  that  the  performance  of  the  NSD  and 
TLS-ESPRIT  estimates  change  with  the  different  values  of  L.  We 
observed  that  the  MSE  starts  decreasing  as  L  increases  for  both 
the  methods  and  for  N=25,  the  best  performance  (min  MSE)  for  the 
NSD  occurs  at  L=15  (s  |n)  and  at  L=12  (s  ^N)  for  TLS-ESPRIT.  The 
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performance  of  the  MNSD  estimates  does  not  seem  to  be  much 
affected  with  the  variation  of  L.  We  report  the  performances  of 
MNSD,  Quinn's  method  (1994)  and  the  best  performance  of  the  NSD 
and  TLS-ESPRIT  along  with  their  CRLB  in  Table  1. 

We  also  performed  a  simulation  study  to  investigate  the 
performance  of  the  different  confidence  intervals  discussed  in 
section  5  with  respect  to  their  average  length  and  coverage 
probabilities.  We  consider  the  simulation  model  (6.1)  with 
white  and  Gaussian  having  error  variance  cr^ .  Results  are 
obtained  for  u  =  0.25Tr,  O.SOtt  and  O.VSrr.  For  each  cj,  simulations 
were  performed  for  <t^  =  0.01,  0.05  and  0.1.  Average  length  of 
the  confidence  intervals  (with  nominal  level  0.90)  and  the 
corresponding  coverage  probabilities  over  100  simulations  are 
reported  for  all  the  methods  in  table  2.  The  bootstrapping 
number  NBOOT  is  taken  as  100  for  both  the  bootstrap  methods. 


7.  CONCLUSIONS 

In  this  article,  we  propose  a  new  non— iterative  method  for 
estimating  the  frequencies  of  the  model  (1.1)  when  the  number  of 
frequencies  is  known  apriori.  If  the  number  of  harmonic 
components  is  unknown,  then  we  can  first  estimate  the  number  of 
frequencies  by  the  method  of  Kundu  (1992)  and  then  use  the 
proposed  method  to  estimate  the  frequencies.  First  we  transform 
the  model  (1.1)  to  an  undamped  superimposed  exponential  signals 
model,  then  use  extended  order  modelling  and  decompose  the  noise 
space  by  singular  value  decomposition  technique.  It  is  further 
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proved  that  the  proposed  non-iterative  technique  yields  estimates 
that  are  strongly  consistent. 

Simulation  results  show  very  satisfactory  performance  of  the 
proposed  method.  The  proposed  method  performs  better  than  the 
best  known  non-iterative  techniques  like  TLS-ESPRIT  (Roy  and 
Kailath;  1989}  and  Quinn's  method  (1994)  even  for  reasonably 
small  sample  sizes  and  at  all  a*  levels.  The  performance  of  MNSD 
almost  attains  the  CRLB  in  the  cases  considered. 

The  choice  of  L  obviously  affects  the  performance  of  the  NSD 
estimates.  Clearly  L  should  be  at  least  M+1,  but  the  natural 
question  is  why  it  should  be  larger  than  that?  Although  no 
theoretical  justifications  have  been  given  in  the  literature,  but 
it  is  observed  that  extended  order  modelling  always  helps  to 
improve  the  performance  of  the  estimators.  Some  heuristic 
justifications  can  be  found  in  Tufts  and  Kumaresan  (1982)  •  It 
seems  more  theoretical  work  is  needed  in  this  direction.  Here  we 
have  observed  that  as  L  increases  the  MSE  starts  decreasing  for 
NSD.  It  reaches  a  minimum  at  L=15  (a  |n)  ,  when  the  sample  size 
is  25.  The  performance  of  the  MNSD  estimates  does  not  seem  to  be 
affected  much  with  variation  in  L. 

Among  the  three  confidence  intervals  for  the  frequencies 
discussed  in  section  5,  the  bootstrap-t  confidence  intervals 
gives  the  highest  coverage  probabilities  although  the  average 
length  of  these  intervals  is  marginally  larger  than  the  other 
two.  It  is  also  obseir/ed  that  the  asymptotic  confidence  interval 
performs  better  than  the  percentile  bootstrap  intervals  in  terms 
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of  shorter  average  length  and  higher  coverage  probabilities  when 
^  ~  O.SOtt.  But  the  percentile  bootstrap  gives  higher  coverage 
probabilities  and  almost  same  length  intervals  when  u)  =  0.2571  or 
0.7577  as  compared  with  the  asymptotic  confidence  intervals.  It 
is  further  observed  that  all  the  three  confidence  intervals  are 
symmetric  about  the  true  parameter  in  the  cases  considered. 
Based  on  the  results  of  the  simulations,  we  recommend  to  use 

bootstrap— t  method  for  finding  confidence  intervals. 

The  research  of  the  authors  is  partially  supported  by  the  Army  Research  Office 
under  the  Grant  DAAH04-96- 1-0082. 
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lauxe  X 


(<J  =  .  25n 

Oi  =  .  SOTT 

(J  =  .  75ti 

CRLB 

0.785398 

3.94548E-06 

1.570796 

3 .41880E-06 

2.356194 

3.03997E-06 

NSD 

0.785401 

4.26988E-06 

1.570874 

4.28355E-06 

2.356232 

4.24039E-06 

MNSD 

0.785388 

3.47565E-06 

1.570840 

3.55599E-06 

2.356200 

3 . 18853E-06 

TLS^ESP 

0.785362 

4.31333E-06 

1.570867 

3.74651E-06 

2.356207 

3.69340E-06 

QUINN 

0.780761 

2.91065E-05 

1.571113 

4.54485E-06 

2.366536 

1. 10624E-04 

CRLB 

0.785398 

7.89097E-06 

1.570796 

6.83761E-06 

2.356194 

6.07997E-06 

NSD 

0.785426 

8.53331E-06 

1.570906 

8 . 57883E-06 

2.356228 

8.50615E-06 

MNSD 

0.785370 

6.93120E-06 

1.570855 

7 . 12682E-06 

2.356217 

6.39456E-06 

TLS-ESP 

0.785353 

8.64388E-06 

1.570896 

7 . 50853E-06 

2.356210 

7.41247E-06 

QUINN 

0.780733 

3.70177E-05 

1.571124 

9.01441E-06 

2.366513 

1.13813E-04 

CRLB 

0.785398 

1.18364E-05 

1.570796 

1.02564E-05 

2.356194 

9.11991E-06 

NSD 

0.785455 

1.27951E-05 

1.570932 

1. 28827E-05 

2.356218 

1.27903E-05 

MNSD 

0.785352 

1. 03701E-05 

1.570865 

1. 07106E-05 

2.356235 

9. 61498E-06 

TLS-ESP 

0.785348 

1.29883E-05 

1.570917 

1. 12857E-05 

2.356211 

1. 11506E-05 

QUINN 


0.780704 

4.49808E-05 


1.571126 

1.34930E-05 


2.366489 

1.17011E-04 


Table  1  (continued) 


(j)  =  .2571 

(J  =  .  5077 

w  =  .  7577 

CRLB 

0.785398 

1.57819E-05 

1.570796 

1. 36752E-05 

2.356194 

1.21598E-05 

NSD 

0.785486 

1.70573E-05 

1.570953 

1.71940E-05 

2.356203 

1.70906E-05 

MNSD 

0.785331 

1.38102E-05 

1.570873 

1.43069E-05 

2.356254 

1.28486E-05 

TLS-ESP 

0.785344 

1.73457E-05 

1.570935 

1. 50783E-05 

2.356212 

1.49056E-05 

QUINN 

0.780674 

5.29854E-05 

1.571122 

1.79825E-05 

2.366466 

1.20221E-04 

CRLB 

0.785398 

1.97274E-05 

1.570796 

1. 70940E-05 

2.356194 

1.51998E-05 

NSD 

0.785519 

2.13213E-05 

1.570972 

2. 15121E-05 

2.356185 

2 . 14056E-05 

MNSD 

0.785313 

1.72385E-05 

1.570880 

1.79152E-05 

2 . 356274 
1.60951E-05 

TLS-ESP 

0.785343 

2.17155E-05 

1.570951 

1.88861E-05 

2.356213 

1.86763E-05 

QUINN 

0.780644 

6.10271E-05 

1.571116 

2.24838E-05 

2.366442 

1.23440E-04 

CRLB 

0.785398 

1.97274E-04 

1.570796 

1.70940E-04 

2.356194 

1.51998E-04 

NSD 

0.787414 

2.24305E-04 

1.571413 

2 . 20710E-04 

2.354701 

2.28583E-04 

MNSD 

0.784367 

1.70489E-04 

1.570889 

1. 92013E-04 

2.357353 

1.75043E-04 

TLS-ESP 

0.785474 

2.31914E-04 

1.571154 

2. 10385E-04 

2.356044 

2.04808E-04 

QUINN 

0.779592 

5.06561E-04 

1.568753 

3.31991E-04 

2.364463 

3.94745E-04 

Table  1  (continued) 


cr2 

(J  =  .257T 

0)  =  .50tt 

0)  =  .  75n 

CRLB 

0.785398 

3.94548E-04 

1.570796 

3.41880E-04 

2.356194 

3.03997E-04 

NSD 

0.790650 

5.39997E-04 

1.571914 

4.52201E-04 

2.351608 

5.40340E-04 

o 

• 

H 

MNSD 

0.783273 

3.47738E-04 

1.570815 

4. 13789E-04 

2.358698 

3.88778E-04 

TLS-ESP 

0.777414 

6.67020E-03 

1.570930 

5.15182E-04 

2.356122 

4.89108E-04 

QUINN 

0.783841 

1.28470E-03 

1.564314 

7.40184E-04 

2.365571 

5.00275E-03 

CRLB 

0.785398 

5.91823E-04 

1.570796 

5. 12821E-04 

2.356194 

4.55995E-04 

NSD 

0.800560 

4.07110E-03 

1.574451 

8. 08582E-04 

2.341074 

4.22644E-03 

1.5 

MNSD 

0.785737 

2.82607E-03 

1.570805 

6.72624E-04 

2.341074 

2.82559E-03 

TLS-ESP 

0.778216 

7.82090E-03 

1.600983 

5. 08289E-02 

2.368486 

1.00645E-02 

QUINN 

0.781995 

2.01329E-03 

1.08511E-02 

1.563569 

1.72613E-03 

8.09113E-03 

2.365616 

8.79353E-03 

2.46323E-02 

Table  2 


Asymptotic 
(Cov.  Prob.) 
Per  -  Boot 
(Cov.  Prob.) 
Bootstrap-t 
(Cov.  Prob.) 


0.005736 

(.82) 

0.006129 

(.88) 

0.006784 

(.89) 


0.005716 

(.85) 

0.005764 

(.84) 

0.006377 

(.84) 


Asymptotic 
(Cov.  Prob.) 
Per  -  Boot 
(Cov.  Prob.) 
Bootstrap-t 
(Cov.  Prob.) 


0.012842 

(.81) 

0.013704 

(.88) 

0.015194 

(.89) 


0.012799 

(.85) 

0.012916 

(.84) 

0.014297 

(.85) 


0.005712 

(.82) 

0.005495 

(.81) 

0.006062 


(.82) 


0.012784 

(.82) 

0.012348 

(.81) 

0.013623 

(.82) 


